function [OF, psi,rhoj,sj,s0j,etaEj,etaDj,vdj, piDprimej, dealer_ids_out, uB, uE, piB, piE, YB, YE, Ecost, uE2, psi_opt, rhoj_opt, W, Wxij, Wvdj, W2] = model_prediction_avg_eta(q, Yt, params, side, estimate, access,vdjSQ,idxj, CF)

global Imax Dmax eta

%sort Yt on id and then dealer 
Yt = sortrows(Yt, [1,16,17]); %Note: this is done in the mainfile as well, should not do anything

%index for dealers
dealer_ids  = Yt(idxj,16);
idxd        = find(ismember(Dmax,dealer_ids));

dealer_ids_out       = NaN*ones(Imax,1);
dealer_ids_out(idxd) = dealer_ids;

if CF>=8
    q=q(idxd);
end 

%label params
sigma   = Yt(1,15);
munu    = params(1);
sigmanu = params(2);
if access==1
c       = params(3);
end 
xij         = Yt(idxj,14);

%distribution functions
f = @(x) normpdf(x,munu,sigmanu);
F = @(x) normcdf(x,munu,sigmanu);
   
%initialize so that all outputes have the same number of dealers
psi        =  NaN*ones(Imax,1);
rhoj       =  NaN*ones(Imax,1);
psi_opt    =  NaN*ones(Imax,1);
rhoj_opt   =  NaN*ones(Imax,1);
s0j        =  NaN*ones(Imax,1);
sj         =  NaN*ones(Imax,1);
etaEj      =  NaN*ones(Imax,1);
etaDj      =  NaN*ones(Imax,1);
piDprimej  =  NaN*ones(Imax,1);
vdj        =  NaN*ones(Imax,1);
uB         =  NaN*ones(Imax,1);
uE         =  NaN*ones(Imax,1);
uE2        =  NaN*ones(Imax,1);
piB        =  NaN*ones(Imax,1);
piE        =  NaN*ones(Imax,1);
YB         =  NaN*ones(Imax,1);
YE         =  NaN*ones(Imax,1);
Ecost      =  NaN*ones(Imax,1);

%use mean theta per dealer (to have a single theta per day)
theta= ones(Imax,1)*mean(Yt(idxj, 5));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%BUYERS  WITH ACCESS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if access==1 && side==1 
    
        %Psi
        psi(idxd,:)     =  sigma*log(sum(exp(1/sigma*(xij+q))));
        psi_opt(idxd,:) =  max(xij+q);
        
        %Platform market share (cond on entry)
        s0j(idxd,:)      = exp(1/sigma*(xij + q)) /(sum(exp(1/sigma*(xij + q))));

        s0jprime(idxd,:) = 1/sigma*(1- s0j(idxd,:)).*s0j(idxd,:);
        
        %Find valueDj
            cut      = @(v)  (psi(idxd,:) + c - eta *(v + xij))./(1-eta) -theta(idxd);
            cutprime = @(v)  s0j(idxd,:)./(1-eta);
            rhoprime = @(v)  -f(cut(v)).*cutprime(v);
            Sprime   = @(v)  sum(f(cut(v)).*cutprime(v)).*s0j(idxd,:) + sum(F(cut(v))).*s0jprime(idxd,:);
            S        = @(v)  sum(F(cut(v))).*s0j(idxd,:);

        if estimate==1
            %Find v s.t. sum(FOC.^2)=0
            funval   = @(v) sum(((1-eta).*(v - (theta(idxd)+cut(v)-xij)).*rhoprime(v) + Sprime(v).*(v - q) - S(v)).^2);
            [vdjout,fval,flag,output] = fmincon(funval, theta(idxd));        
            vdj(idxd,:) = vdjout;
        
        elseif estimate==0         
            vdj(idxd,:)     =vdjSQ(:,idxd)';%use the estimated value of the dealer
        end 

        %Bilateral market share 
        rhoj(idxd,:)     = 1-F(cut(vdj(idxd,:)));
        rhoj_opt(idxd,:) = 1-F((psi_opt(idxd,:) + c - eta *(vdj(idxd,:) + xij))./(1-eta) -theta(idxd));

        %Platform market share unconditional on venue choice
        sj(idxd,:)      = sum((1-rhoj(idxd,:))).* s0j(idxd,:);

        if CF==8 || CF==9 || CF==10 || CF==11
            s0j(idxd,:)=0;
            s0j(max(xij+q)==xij+q,:)=1;
            sj(idxd,:)= s0j(idxd,:).* sum(F(psi(idxd,:) +c - theta(idxd)));
        end 
        
        %Elasticities 
        etaEj(idxd,:) = Sprime(vdj(idxd,:))./S(vdj(idxd,:)).*q ;
        etaDj(idxd,:) = rhoprime(vdj(idxd,:))./rhoj(idxd,:).*q;
        
        %Marginal profit in bilateral
        piDprimej(idxd,:) = (1-eta).*(vdj(idxd,:) - (theta(idxd)+cut(vdj(idxd,:))-xij)).*rhoprime(vdj(idxd,:));


        %OBJECTIVE FUNCTION 
         OF0  =  sum(((1-eta).*(vdj(idxd,:)  - (theta(idxd)+cut(vdj(idxd,:) )-xij)).*rhoprime(vdj(idxd,:)) + Sprime(vdj(idxd,:)).*(vdj(idxd,:) - q) - S(vdj(idxd,:))).^2);
        if isnan(OF0)==0 
            OF = OF0;
        else 
            OF = 10^100;
        end 

 
        %EXPECTED PROFIT OF INVESTOR      
        %a) from bialteral trade
        uB(idxd,:) = eta.*(vdj(idxd,:) -(theta(idxd)+munu-xij));
        
        %b) from platform trade 
        a             = cut(vdj(idxd,:));
        Mills         = -f(a)./F(a);  
        ExpfunT       = munu + sigmanu.*Mills ;
        uE(idxd,:)    = -ExpfunT + (1-rhoj(idxd,:)).*(-theta(idxd)+ (psi(idxd,:)) );
        uE2(idxd,:)   = -ExpfunT + (1-rhoj(idxd,:)).*(-theta(idxd) + sum((s0j(idxd,:).*(q+xij)))); %exclude sigma*epsilon
        Ecost(idxd,:) = (1-rhoj(idxd,:)).*c;


        %EXPECTED YIELD
        if F(a)~=1     %some trade bilateral
            YB(idxd,:) = rhoj(idxd,:).*(eta.*vdj(idxd,:)+(1-eta).*(theta(idxd)-xij)) + (1-eta).*(munu + sigmanu.*(f(a)./(1-F(a))));
        elseif F(a)==1 %all trade on the platform
            YB(idxd,:)=0;
        end    
        YE(idxd,:) = (1-rhoj(idxd,:)).* sum((s0j(idxd,:).*q)); 
        

        %EXPECTED PROFIT FOR THE DEALER
        %Profit from bilateral trades 
        piB(idxd,:) = (1-eta).*(vdj(idxd,:) -theta(idxd) +xij).*(1-F(a)) ...
                    - (1-eta).*(munu + sigmanu*(f(a)./(1-F(a))));
                
        %Profit from platform trades
        piE(idxd,:) = sj(idxd,:).*(vdj(idxd,:)-q);
        
        
        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)+xij-(theta(idxd)+munu)));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)));
   
        %Welfare when shrinking dealer's value towards theta
        distance = vdj-theta;
        vdj2     = theta+(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));

         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%SELLERS  WITH ACCESS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif side==-1 && access==1
            
        %psi
        psi(idxd,:)     = -sigma*log(sum(exp(1/sigma*(xij-q))));
        psi_opt(idxd,:) = -max(-q+xij);
        
        %Platform market share (cond on entry)
        s0j(idxd,:)      = exp(+ 1/sigma*(xij - q)) /(sum(exp(+1/sigma*(xij - q))));
        s0jprime(idxd,:) = -1/sigma*(1- s0j(idxd,:)).*s0j(idxd,:);
        
        %Find valueDj
            cut      = @(v)  (psi(idxd,:) - c + eta *(-v + xij))./(1-eta) -theta(idxd);
            cutprime = @(v)  s0j(idxd,:)./(1-eta);
            rhoprime = @(v)  f(cut(v)).*cutprime(v);
            Sprime   = @(v)  -sum(f(cut(v)).*cutprime(v)).*s0j(idxd,:) + sum((1-F(cut(v)))).*s0jprime(idxd,:);
            S        = @(v)  sum((1-F(cut(v)))).*s0j(idxd,:);

       if estimate==1
       %Find v s.t. sum(FOC.^2)=0
            funval   = @(v) sum(((1-eta).*((theta(idxd)+cut(v)+xij) - v).*rhoprime(v) + Sprime(v).*(q-v) + S(v)).^2);
            [vdjout,fval,flag,output] = fmincon(funval, theta(idxd));        
            vdj(idxd,:) = vdjout;
        
        elseif estimate==0         
            vdj(idxd,:)     =vdjSQ(:,idxd)';%use the estimated value of the dealer
       end 

        %Bilateral market share 
        rhoj(idxd,:)     = F(cut(vdj(idxd,:)));
        rhoj_opt(idxd,:) = F((psi_opt(idxd,:) - c + eta *(vdj(idxd,:) + xij))./(1-eta) -theta(idxd));

        %Platform market share unconditional on venue choice
        sj(idxd,:)      = sum((1-rhoj(idxd,:))).* s0j(idxd,:);

        %Elasticities 
        etaEj(idxd,:) = Sprime(vdj(idxd,:))./S(vdj(idxd,:)).*q ;
        etaDj(idxd,:) = rhoprime(vdj(idxd,:))./rhoj(idxd,:).*q;
      
        %Marginal profit in bilateral
        piDprimej(idxd,:) = (1-eta).*( (theta(idxd)+cut(vdj(idxd,:))+xij)-vdj(idxd,:)).*rhoprime(vdj(idxd,:));

        %OBJECTIVE FUNCTION 
         OF0  =  sum(((1-eta).*((theta(idxd)+cut(vdj(idxd,:))+xij)-vdj(idxd,:)).*rhoprime(vdj(idxd,:)) + Sprime(vdj(idxd,:)).*(q-vdj(idxd,:)) + S(vdj(idxd,:))).^2);
        if isnan(OF0)==0 
            OF = OF0;
        else 
            OF = 10^100;
        end 

        %EXPECTED PROFIT OF INVESTOR      
        %a) from bialteral trade
        uB(idxd,:) = eta.*((theta(idxd)+munu+xij)-vdj(idxd,:));
        
        %b) from platform trade 
        a             = cut(vdj(idxd,:));
        Mills         = f(a)./F(a);  
        ExpfunT       = munu + sigmanu.*Mills ;
        uE(idxd,:)    = ExpfunT + (1-rhoj(idxd,:)).*(theta(idxd)- (psi(idxd,:))); %include sigma*epsilon
        uE2(idxd,:)   = ExpfunT + (1-rhoj(idxd,:)).*(theta(idxd)- sum((s0j(idxd,:).*(q-xij)))); %exclude sigma*epsilon
        Ecost(idxd,:) = (1-rhoj(idxd,:)).*c;
   
        %EXPECTED YIELD 
        if F(a)~=0 %some trades billaterally
            YB(idxd,:) = rhoj(idxd,:).*(eta.*vdj(idxd,:) +(1-eta).*(theta(idxd) +xij)) + (1-eta).*(munu - sigmanu.*(f(a)./F(a)));
        elseif F(a)==0 %no one trades billaterally
            YB(idxd,:)=0;
        end 
        YE(idxd,:) = (1-rhoj(idxd,:)).* sum((s0j(idxd,:).*q)); 
         
        
        %EXPECTED PROFIT FOR THE DEALER 
        %Profit from bilateral trades 
        piB(idxd,:) = (1-eta).*(theta(idxd) +xij - vdj(idxd,:)).*F(a) ...
                    + (1-eta).*(munu - sigmanu*(f(a)./F(a)));
                
        %Profit from platform trades
        piE(idxd,:) = sj(idxd,:).*(q-vdj(idxd,:));
              

        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*((theta(idxd)+munu)-vdj(idxd,:)+xij));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(-vdj(idxd,:)));
              
        distance = vdj-theta;
        vdj2     = theta+(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
       

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%BUYERS  WITHOUT ACCESS       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif access==0 && side==1  
    
        %Use the estimated value of the dealer
        vdj(idxd,:)     =vdjSQ(:,idxd)';
    
        %PROBABILITIES
        rhoj(idxd,:) = ones(size(idxd,2),1); 
        sj(idxd,:)   = zeros(size(idxd,2),1);  
        s0j(idxd,:)  = zeros(size(idxd,2),1);  
       
        %EXPECTED PROFIT OF INVESTOR    
        uB(idxd,:) = eta.*(vdj(idxd,:) -(theta(idxd)+munu-xij));
        uE(idxd,:) =0; 
        Ecost(idxd,:)=0; 
        YE(idxd,:) =0; 
        YB(idxd,:) = eta.*vdj(idxd,:)+(1-eta).*(theta(idxd)-xij) ; 
 
        %EXPECTED PROFIT FOR THE DEALER
        %Profit from bilateral trades 
        piB(idxd,:) = (1-eta).*(vdj(idxd,:) -theta(idxd) -munu +xij);
        %Profit from platform trades
        piE(idxd,:) = 0.*vdj(idxd,:);
        OF=NaN;
        
        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)+xij-(theta(idxd)+munu)));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)));

        %Shrinking dealer's value towards theta
        distance = -vdj+theta;
        vdj2     = theta-(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
         

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%SELLERS  WITHOUT ACCESS   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif side==-1 && access==0  
    
        %Use the estimated value of the dealer
        vdj(idxd,:)     =vdjSQ(:,idxd)';
        
        %Probability to trade OTC
        rhoj(idxd,:)= ones(size(idxd,2),1);  
        sj(idxd,:)  = zeros(size(idxd,2),1);  
        s0j(idxd,:) = zeros(size(idxd,2),1);  
        
        uB(idxd,:) = eta.*((theta(idxd)+munu+xij)-vdj(idxd,:)) ; 
        uE(idxd,:) =0; 
        Ecost(idxd,:) =0; 
        YE(idxd,:) =0; 
        YB(idxd,:) = (1-eta).* vdj(idxd,:)  + eta.*( theta(idxd) + munu + xij);
        
        %Profit from bilateral trades 
        piB(idxd,:) = (1-eta).*(theta(idxd) + munu +xij  - vdj(idxd,:));
        %Profit from platform trades
        piE(idxd,:) = 0.*vdj(idxd,:);
             
        %OF set to NaN (never maximized over!)
        OF=NaN;

        %Welfare
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*((theta(idxd)+munu)-vdj(idxd,:)+xij));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(-vdj(idxd,:)));

        %Shrinking dealer's value towards theta
        distance = -vdj+theta;
        vdj2     = theta-(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
     
end
